C************************************************************* C Test program for subroutine LSEMA by Stephen Kirkup C************************************************************* C C Version 2. Revised test problems. Changes C in interface with LSEM3. C School Engineering, University of Central Lancashire C www.uclan.ac.uk C smkirkup@uclan.ac.uk C http://www.researchgate.net/profile/Stephen_Kirkup C C This open source code can be found at C www.boundary-element-method.com/fortran/LSEM3_T.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C C*************************************************************** C This program is a test for the subroutine LSEMA. The program computes C the solution to the Laplace equation C C __ 2 C \/ {\phi} = 0 C C exterior to a (set of) thin shell(s). C C Subroutine LSEMA accepts a description of the shell of the domain C and the position of the exterior points where the solution ({\phi}) C is sought, the shell condition and returns the solution ({\phi} C and v) on S and the value of {\phi} at the exterior points. C C C The test problems C ----------------- C C In this test the domain is the exterior to two circular plates of C radius 1 (metre) and 0.1 apart. The solution to the problem with a C Dirichlet shell condition ({a}=1, {b}=0) with C {\phi}=f=1.0 on the upper plate and {\phi}=f=0.0 on the lower plate. C C We would expect a simple gradient in potential between the two C plates and this is investigated in this test program. C C The plates are described by a set of NH=20 panels of equal length C along the generator. The shell solution points are the centres of the C generator of the panels. C The solution is sought the points (0,0,0.025), (0,0,0.05), C (0,0,0.075). C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNH : The limit on the number of shell panels. C integer MAXNV : The limit on the number of vertices. C integer MAXNPE : The limit on the number of exterior points. C External modules related to the package C --------------------------------------- C subroutine LSEMA: Subroutine for solving the exterior Laplace C equation. (file LSEMA.FOR contains LSEMA and subordinate routines) C subroutine L3ALC: Returns the individual discrete Laplace integral C operators. (file L3ALC.FOR contains L3ALC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS2.FOR ) C file GEOM2D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LSEMAT IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of panels INTEGER MAXNH PARAMETER (MAXNH=32) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=MAXNH) C Limit on the number of points exterior to the shell, where C the solutions are sought INTEGER MAXNPE PARAMETER (MAXNPE=6) C Geometrical description of the shell(s) C Number of panels and counter INTEGER NH,IH C Number of central points (on H) and counter INTEGER NHP,IHP C Number of vetices and counter INTEGER NV,IV C Index of nodal coordinate for defining shell(s) standard unit is C metres) REAL*8 VERTEX(MAXNV,2) C The two nodes that define each panel on the boundaries INTEGER HELV(MAXNH,2) C The points exterior to the shell(s) where the solutions C are sought and the directional vectors at those points. C [Only necessary if an exterior solution is sought.] C Number of exterior points and counter INTEGER NPE,IPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,2) C Data structures that are used to define each test problem in turn C and are input parameters to LSEMA. C HA(j) is assigned the value of {a} at the centre of the C j-th panel. REAL*8 HA(MAXNH) C HB(j) is assigned the value of {b} at the centre of the C j-th panel. REAL*8 HB(MAXNH) C HF(j) is assigned the value of f at the centre of the j-th panel. REAL*8 HF(MAXNH) C HA(j) is assigned the value of A at the centre of the C j-th panel. REAL*8 HAA(MAXNH) C HB(j) is assigned the value of B at the centre of the C j-th panel. REAL*8 HBB(MAXNH) C HF(j) is assigned the value of F at the centre of the j-th panel. REAL*8 HFF(MAXNH) C The incident field C The incident potential on the shell(s) REAL*8 HIPHI(MAXNH) C The derivative of the incident potential on the shell(s) REAL*8 HIVEL(MAXNH) C The incident potential at the exterior points REAL*8 PEIPHI(MAXNH) C Validation and control parameters for LSEMA C Switch for particular solution LOGICAL LSOL C Validation switch LOGICAL LVALID C The maximum absolute error in the parameters that describe the C geometry of the shell. REAL*8 EGEOM C Output from subroutine LSEMA C The average potential at the centres of the panels REAL*8 PHIAV(MAXNH) C The difference in potential at the centres of the panels REAL*8 PHIDIF(MAXNH) C The average derivative of the potential at the centres of the C panels REAL*8 VELAV(MAXNH) C The difference in the derivative of the potential at the centres C of the panels REAL*8 VELDIF(MAXNH) REAL*8 PEPHI(MAXNPE) C Working space REAL*8 AMAT(2*MAXNH,2*MAXNH) REAL*8 BMAT(2*MAXNH,2*MAXNH) REAL*8 L_EH(MAXNPE,MAXNH) REAL*8 M_EH(MAXNPE,MAXNH) C Further information from system solution GLS INTEGER*4 PERM(2*MAXNH) LOGICAL XORY(2*MAXNH) REAL*8 C(2*MAXNH) C Work Space REAL*8 WKSPC1(2*MAXNH) REAL*8 WKSPC2(2*MAXNH) REAL*8 WKSPC3(2*MAXNH) C Counter through the x,y coordinates INTEGER ICOORD C INITIALISATION C -------------- C Describe the nodes and panels that make up the shell C :The plates are circles, each divided into 10 panels. VERTEX and C : HELV are defined out from the centres of the plates so that the C : normal to the plates is assumed to be upward (+z direction). C :Set up nodes C : Set NH, the number of panels NH=20 C : Set NV, the number of vertices (equal to the number of panels) NV=22 C : Set coordinates of the nodes DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,22) * / 0.000D0, 0.000D0, * 0.100D0, 0.000D0, * 0.200D0, 0.000D0, * 0.300D0, 0.000D0, * 0.400D0, 0.000D0, * 0.500D0, 0.000D0, * 0.600D0, 0.000D0, * 0.700D0, 0.000D0, * 0.800D0, 0.000D0, * 0.900D0, 0.000D0, * 1.000D0, 0.000D0, * 0.000D0, 0.100D0, * 0.100D0, 0.100D0, * 0.200D0, 0.100D0, * 0.300D0, 0.100D0, * 0.400D0, 0.100D0, * 0.500D0, 0.100D0, * 0.600D0, 0.100D0, * 0.700D0, 0.100D0, * 0.800D0, 0.100D0, * 0.900D0, 0.100D0, * 1.000D0, 0.100D0/ C :Describe the panels that make up the two shells C : Set NH, the number of panels NH=20 C : Set nodal indices that describe the panels of the shells. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in HELV dictates that the normal is outward from the C : shell into the domain. DATA ((HELV(IH,ICOORD),ICOORD=1,2),IH=1,20) * / 1, 2, * 2, 3, * 3, 4, * 4, 5, * 5, 6, * 6, 7, * 7, 8, * 8, 9, * 9, 10, * 10, 11, * 12, 13, * 13, 14, * 14, 15, * 15, 16, * 16, 17, * 17, 18, * 18, 19, * 19, 20, * 20, 21, * 21, 22 / C Set the points in the domain where the solutions are sought, PEXT. NPE=3 DATA ((PEXT(IPE,ICOORD),ICOORD=1,2),IPE=1,3) * / 0.000D0, 0.025D0, * 0.000D0, 0.050D0, * 0.000D0, 0.075D0/ C The number of points on the shell is equal to the number of C panels NHP=NH C :In the test problem the plate 1 has C {\phi}= 1 C : and plate 2 has C {\phi}= 0 DO 250 IHP=1,NHP/2 HA(IHP)=1.0D0 HB(IHP)=0.0D0 HF(IHP)= 0.0D0 HAA(IHP)=1.0D0 HBB(IHP)=0.0D0 HFF(IHP)= 0.0D0 250 CONTINUE DO 260 IHP=NHP/2+1,NHP HA(IHP)=1.0D0 HB(IHP)=0.0D0 HF(IHP)= 0.0D0 HAA(IHP)=1.0D0 HBB(IHP)=0.0D0 HFF(IHP)= 1.0D0 260 CONTINUE C There is no incident field DO 640 IHP=1,NHP HIPHI(IHP)=0.0D0 HIVEL(IHP)=0.0D0 640 CONTINUE DO 650 IPE=1,NPE PEIPHI(IPE)=0.0D0 650 CONTINUE C Set up validation and control parameters C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LSEMA LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Set up particular a,b,f,A,B,F functions for this wavenumber C and type of shell condition CALL LSEMA(MAXNV,NV,VERTEX,MAXNH,NH,HELV, * MAXNPE,NPE,PEXT, * HA,HB,HF,HAA,HBB,HFF, * HIPHI,HIVEL,PEIPHI, * LSOL,LVALID,EGEOM, * PHIDIF,PHIAV,VELDIF,VELAV,PEPHI, * AMAT,BMAT,L_EH,M_EH, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='LSEMA.OUT') WRITE(20,*) 'LSEMA: Computed solutions in the domain' WRITE(20,*) WRITE(20,*) ' P(1) P(2) Computed Exact/App' DO 800 IPE=1,NPE WRITE(20,999) PEXT(IPE,1),PEXT(IPE,2),PEPHI(IPE), * PEXT(IPE,2)*10.0D0 800 CONTINUE 999 FORMAT(4F10.4) CLOSE(20) END